Nonlinear EHD instability of two viscoelastic fluids under the influence of mass and heat transfer

This study attempts to provide an approach to studying the nonlinear stability of a vertical cylindrical interface between two Oldroyd-B prototypes. An unchanged axial electric field influences the system, and porous medium, and the effects of heat and mass transfer (MHT) are considered. Hsieh's modulation and the viscous potential flow (VPT) are used to abbreviate the mathematical analysis. The viscoelastic Oldroyd-B model significant role in geothermal, engineering and industrial enhancement motivated us to carry out this in-depth investigation. The methodology of the nonlinear technique depends mainly on solving the linear equations of motion and applying the appropriate nonlinear boundary conditions. Numerous non-dimensional physical numbers are exposed using a non-dimensional technique. The stability conditions are theoretically achieved and numerically verified. As a limiting case, the linear dispersion equation is accomplished, and a set of stability diagrams is reachable. Together with the nonlinear stability method, a Ginzburg–Landau equation is derived. Subsequently, both theoretical and numerical methods are used to achieve the nonlinear stability criteria. Furthermore, a precise perturbed approach for surface deflection is achieved theoretically and numerically using the Homotopy perturbation method and the extended frequency conception. Along with the linear approach, it is found that the structure becomes unstable by the Laplace, Reynolds, Weber, and elasticity quantities as well as the linear MHT parameter. Furthermore, the stability zones are enhanced in the nonlinear instability approach.

Electrohydrodynamics (EHD) is an area of Fluid Mechanics regarding electrical strength impacts 1 . It could be regarded as an issue that deals with moving fluids in an electric field. Fundamentally, it is a combination of these two areas because many of the most exciting challenges in EHD include both the influence of liquid motion and the impact of electric strength. On the other side, EHD incorporates a complicated interaction of internal, viscous, and electric forces. It produces attractive phenomena in processing equipment in inkjet printing, drops sparing, and other applications. Rayleigh was the first to study charged droplet instability over a century ago; therefore, EHD has been a hot topic. A review article was constructed to give an integrated picture of EHD 2 . The interface stability separating two layers of dielectric liquids in the existence of an applied voltage in addition to imposed fluid fields was examined. Several pioneering works were done, and many researchers have been studying the EHD stability of the interface dividing two superimposed liquids in both bounded and unbounded channels 3 . Using the surface-connected prototype, the EHD instability of two superimposed liquids in a channel was scrutinized 4 . The EHD stability of an incompressible viscous fluid jet was presumed by a uniform longitudinal electric strength, using a weakly nonlinear approach 5 . The nonlinear boundary conditions were employed along with the linear fundamental equations of motion. The stability of an interface separating different liquids with electrical relaxation time was investigated in depth. The nonlinear EHD of capillary surface waves was investigated 6 . The effect of a longitudinal periodic field on a streaming flow passing through three coaxial vertical cylinders was investigated 7 . Coupled Mathieu equations were conducted in their approach. They considered a standard example of surface displacement. The behavior of an electric strength on two cylindrical boundaries was studied 8 . Through complex quantities, the method produced differential equations with damped terms. For more convenience, these equations were combined in light of the symmetrical and anti-symmetric perturbations. Three coupled magnetic fluid interfaces were analyzed for nonlinear azimuthal stability 9 . Their nonlinear method was constructed primarily on resolving the linear regulating equations of motion and applying the associated nonlinear boundary conditions. A streaming dielectric liquid jet EHD stability was conducted. The system was examined to see if a constant axial electrostatic field had entered it 10 . From their process, several non-dimensional numbers were produced. Their research gave researchers a strong foundation for understanding how viscous liquid jets break up and become unstable in the presence of electric strength. Because of how much the electric strength influences things, the current technique is employed in this direction.
Today, there is a growing and substantial quantity of literature on Newtonian fluid cylinder instability. Viscoelastic liquids are non-Newtonian liquids that demonstrate both viscous and elastic properties under certain conditions. Walters, Rivlin-Ericksen, Maxwell, Kelvin, Oldroyd, and others suggested many forms of viscoelastic models. Fluid flows in geophysics, such as polymer distribution, particle filtration from fluids, and lubrication, are frequently represented as viscoelastic flows. The complicated structure of the fundamental performance of such fluids makes further complicated things for viscoelastic jets 11 . The non-Newtonian liquid takes part in a wide range of productions and scientific disciplines. The Oldroyd-B is considered a good model to investigate further significant experiments among several prototypes that have been employed to illustrate the performance of non-Newtonian liquids. This model was the most successful for comparing the reactions of numerous viscoelastic liquids, since it can predict stress relaxation, creep, and normal stress changes. It covers the Maxwell and viscous liquid simulations as specific cases. Industries are interested in the motion of fluids in the surroundings of a moving body. One of the really significant and fascinating movements was the flowing through revolving cylinders which have different uses in the food sector. Non-Newtonian liquids possess mechanics that are both practically and theoretically significant. The subject was theoretically interesting since the linear instability evaluation of these liquids was helpful in highlighting the properties under certain conditions 12 . When HMT plays a significant role in establishing stream stability, injector-sparing qualities were examined. The Kelvin-Helmholtz Scientific Reports | (2023) 13:357 | https://doi.org/10.1038/s41598-023-27410-z www.nature.com/scientificreports/ attachments breaking apart have been found in prior research. These processes were related to the development of disturbances waves on the liquid-air interface, which occurs downstream of the injector and, once a critical waveform was achieved, cause the liquid jet to fragment into ligaments and/or droplets 35 and 36 . The major objective of the current study is to help readers to look for appropriate responses to the following questions: 1. What are the criteria of the linear stability approach? 2. How many physical non-dimensional numbers are present throughout the linear technique? 3. What are the effects of the nondimensional physical numbers on the stability configuration? 4. What is the approach of the nonlinear stability sense? 5. What is the formula for the approximate solution to the interface displacement?
For ease of understanding, the outstanding physical configuration will be established as follows: Section "Methodology of the model" introduces the approach of the physical prototype. The fundamental controlling equations of motion and the relevant nonlinear border restrictions are involved in this Section. In section "Linear stability methodology", the relevant linear stability methodology is introduced. Additionally, a set of graphs is drawn in order to depict the influences of the physical non-dimensional factors in the instability information. The nonlinear aspects in addition to the relevant numerical estimation are introduced in section "Nonlinear stability methodology". An approximate form of the interface displacement via the HPM and the expanded nonlinear frequency approach is provided in section "Estimated surface departure formulation". Furthermore, two dimensions of the interface displacement are displayed. Finally, in section "Concluding remarks", concluding remarks are formed based on the main ideas.

Methodology of the model
Two homogeneous, incompressible, and dielectric-flowing liquids make up the structure under consideration. The theoretical system consists of three cylindrical surfaces. The inner cylinder is a solid one, and it is raised to a uniform temperature. The flow is governed by Darcy's law throughout porous mediums. For ease of reference, the permeability of the two fluids is simply referred to as. Scientific evidence supports Darcy's law, an experimental hypothesis that describes the creeping flow of Newtonian liquids in porous media. A non-Newtonian viscoelastic fluid obeying an Oldroyd-B flow occupies the inner cylindrical annular liquid. This liquid is characterized by certain constants, which are addressed as retardation and relaxation times in addition to the viscosity parameter. Viscoelasticity is vital for understanding how blood flow behaves, and the constitutive model of Oldroyd-B fluid can be applied in this area. It contains relaxation and retardation times, as two important characteristics of viscoelastic fluid. Restoring an unstable system to equilibrium is typically referred to as relaxing. A relaxation time can be used to categorize each relaxation step. Electrical conductivity and the dielectric relaxation time are strongly connected. It is an indicator of how long it takes for an electrical charge to be neutralized in a semiconductor. Metals have a short relaxation time, whereas semiconductors and insulators can have a long one. Retardation, also known as "delay of the elasticity," is the delayed response to an applied force or stress. Ideal elastic materials exhibit an instantaneous deformation with the application of a tension resembling a leap and an instantaneous reformation upon the removal of the load. For instances of viscoelastic materials, this elastic behavior happens after a specific amount of time. Numerous polymer industries, blood flow, etc. use relaxation/retardation time as a crucial characteristic of viscoelastic fluid flows. A portion of the energy in the cardiovascular system is stored because the blood is elastic, a portion is lost as heat because the blood is viscous, and a portion is involved in blood mobility 37 . A viscous gas occupies the outer cylinder and possesses another viscosity constant. In the undisturbed situation, these two phases are divided by a hypothetical cylindrical interface. This interface is assumed to have another uniform temperature. Finally, the inner cylinder is a solid one and it is raised to distinct constant temperatures. The two phases basically have uniformly various velocities. A longitudinal unchanged electric strength is portrayed in the system. Therefore, the cylindrical coordinates are more convenient for this clarification. The surface tension parameter is taken into the explanation. To facilitate the mathematical manipulation, the MHT is analyzed in accordance with Hsieh's explanation 38 and 39 . This oversimplification results in a single parameter in the linear stability approach, as will be seen later. The nonlinear methodology, on the other side, yields only two parameters to control MHT. Accordingly, this approach would substantially facilitate the relationship experiment. Alongside the negative z− path, the gravitational force is considered. The fluid jet is stable for all asymmetric types, but unstable for the axisymmetric way, as demonstrated by several researchers such as Chandrasekhar 40 . Consequently, the axisymmetric mode is the most interesting one in this study. Therefore, the radiating symmetry is merely reflected. Henceforth, all fundamental amounts will be unrestricted with the angular position. A similar problem with the Walters-B fluids has been recently analyzed 41 . For more explanation, Fig. 1 describes the physical and theoretical model.
The distributed cylindrical surface formulation can be expressed as: where η(z; t) is the surface movement. Consequently, the interface may be given following a small deviation from the equilibrium situation. It could be characterized as: It is possible to formulate the unit normal vector as: where P is the hydrostatic pressure, I denotes the identified tensor, and τ symbolizes an additional tensor that can be described as: The strain tensor rate is also identified as: where N = ∇v.
The Oldroyd-B model is provided by incorporating the traditional Newtonian prototype as a particular situation 1 = 2 = 0 . Furthermore, it is shortened to the Maxwell model at 2 = 0 . Darcy's law states that the decreased pressure is caused by a frictional drag in a viscous Newtonian liquid flowing at a lower speed and directly proportional to velocity. The relationship between pressure drops and the speed of a viscoelastic liquid in a permeable medium may be characterized as previously given 44 in contrast to the Oldroyd-B constitutive relationships.
The localized conservation of momentum balancing is provided by the evidence of the balancing of forces working on a volume amount of liquid.
Because the pressure gradient in Eq. (7) indicates that the flow resistance in the volume of the permeable substance can be measured, the movement resistance provided by the solid matrix is indicated by ϒ . The vector ϒ can be determined from Eq. (7) to fulfill the following equation as 45 : (4) τ vis = −PI + τ ,   (11) can now be rewritten as The Darcy coefficient is equal to ν l = µ l /ξ. The general stress tensor is constructed of two elements. The earliest element is the viscoelastic tensor of the Oldroyd-B prototype which is produced in Eqs. (4) and (5). The mixture of Eqs. (4) and (5) generates: The continuity equation becomes Based on the basic assumptions of the VPF, one might presume that the liquids are irrotational. Accordingly, the velocity potentials exist as the distribution φ j (r, z; t) such that Owing to the incompressibility circumstance, the potential φ j should verify the following Laplace equation: Many investigators have proven that the solutions of distribution functions are typically based on an understanding of natural conditions 9 . As a result of the standard mode approach, one may presume The solution of Eq. (16) then becomes where A j (t) and B j (t) are unspecified time-dependent formulas. They can be calculated employing nonlinear border circumstances. The modified Bessel functions of zero order of the first and second kind, correspondingly, are signified by I 0 (kr) and K 0 (kr).
The pressure is calculated by explaining Eq. (8). Following the approach of Bernoulli's equation, one can find the principal momentum equation by direct integration.
The established Maxwell's formulae must be characterized owing to the overall electrostatic impact on the problem at hand; for instance, see Melcher 47 . Melcher purported a comprehensive book that contained a completed inspection of the surface-coupled EHD and Magnetohydrodynamics (MHD) structures. In the contemporary case, just the impression of a tangential electrostatic strength is reflected. Consequently, in the present study, the presence of magnetic strength can be ignored. Depending on all these assumptions, along with the quasi-static approximation, Maxwell's formulae are reduced to: and Subsequently, the electric strength E j could be denoted with a scalar function ψ j (r, z; t) , as  (20) and (21), the potential of the electric field follows Laplace's equation: As stated in Eq. (17), the electric potentials could be expressed as: where C j (t) and D j (t) are unspecified time-dependent effects that will be assessed later.
Nonlinear border criteria. At the solid cylinders (r = a and r = b).
• The fluid normal velocities should always be eliminated; hence, the following prerequisites should be met.: and • The tangential elements of the electric potential have to be removed, which means and At the free cylindrical boundary r = R + η(z, t).
• The preservation of mass and energy through the surface of separation was achieved by using the abbreviated modulation, which was first suggested by Hsieh 38 and 39 . This created the following requirements.: where � * � = * 2 − * 1 indicates, correspondingly, the difference between the interior and exterior liquid phases.
• According to Hsieh 38 and 39 ., it is expected that the immediate location of the interface has a significant impact on how much latent heat is released. Consequently, the interface requirement for heat exchange is expressed as: where f (η) signifies the net heat flow from the interface and may be formulated as: When expanding f (η) in the around η = 0 , one gets The quantity f (0) denotes the net heat flow from the interface into the liquid zones. Since f (0) is assumed to be zero, the following relationship can be derived: where G is constant, signifying that the heat flux throughout the interface connects the two liquids. It was considered to be identical in the equilibrium condition.
Equations (1), (31), and (32) are now substituted into Eq. (30), yielding the following equation: . and • At the border dividing the two cylindrical liquids, the change in the tangential elements of the electric strength is continuous. Therefore, one finds • At the connection, the continuously normal electric fluctuation provides Replacing Eqs. (18) and (24) into Eqs. (25)- (29) and (35)- (36), it follows that the solutions φ and ψ might be exemplified in the form and where herein I 1 (kR) and K 1 (kR) are the modified Bessel differential equations of first and second kind, respectively. The total stress tensor 48 is provided as follows: for the liquid phase, one gets where τ rr = 2µ l Mφ lrr , τ rz = 2µ l Mφ lrz , τ zz = 2µ l Mφ lzz , and M = µ l 1 + 1 for the gas phase, one finds σ elc ij for the two phases yields σ vis ij = −P g I. www.nature.com/scientificreports/ • The interface tension-induced discontinuities in the normal stress tensor are the source of the residual border criterion. It could be represented as follows.: where F is the complete force on the interface that is described as: where n r and n z are the elements of the unit normal vector n and T is the quantity of the interfacial tension. Therefore, the nonlinear characteristic equation may be expressed as follows: where G(η) represents all the nonlinear terms in the displacement function η . All the constants are supplied in the Appendix to follow the content straightforwardly. It should be emphasized that left-hand side of Eq. (47) demonstrates the linear dispersion relation, allowing us to explain the motion of the cylindrical interface of the linear case. The discussion of this approach is hence the focus of the next Section.
Abbreviation. In the subsequent table, the characters L and signify the internal and external cylindrical fluids, respectively. where m,l,t and K denote mass, length, time, and Kelvin, respectively

Linear stability methodology
It is a respectable impression to look at the stability profile from a linear perspective before moving to the general case. When the nonlinear elements are eliminated from Eq. (47), the dispersion equation is provided for the linear case. The linear dispersion connection can therefore be expressed as: A consistent wave train solution to Eq. (48) can be presumed with the standard mode approach. Therefore, we may write where γ signifies a wave train little amount.
When γ has a nontrivial solution in Eq. (48), the dispersion relation grows to be nontrivial. where A linear dispersion relation that distributes along the two cylindrical flowing fluids along with the fundamental axis is represented by Eq. (50). It contains of the wave number and growth rate that correspond to all of the problem physical constraints. According to the linear stability theory, the growth rate behavior determines whether a period of time is stable or unstable. If the imaginary part is positive, the disruption will spread with time. Consequently, the system progresses linear instability. The scheme would become linearly stable when the imaginary portion of is negative since it will cause a reduction in the disturbance. In these situations, the stability or instability depends on how the imaginary part of the system behaves and whether it oscillates for any time. It should be noted that the real part of the frequency of the surface wave has no implication in the stability configuration. This section major focus on examining the system stability under a linear approach. Furthermore, the nonlinear stability analysis has different criteria rather than the linear one. Simply, the linear modes analysis does not valid throughout the linear sense. At this point, the Routh-Hurwitz theory is applied to achieve the instability criteria of the dispersion relationship (50), see Peña 49 . Subsequently, the following are the stability requirements: and The analysis shows that β 1 is independent of the electric strength E 2 0 . Typically, the stability profile needs to display E 2 0 against the wave numeral. To verify the implication of condition (51), all the character values should satisfy this condition, in order that this condition can become automatically satisfied. On the other hand, the condition (52) can be arranged to be formulated in the following form: , the coefficients b 21 and b 22 are known from the context.
Before performing the numerical calculation, the stability requirement given by inequalities (51) and (53) can be represented in satisfactory formulae. To this end, they must be written in a dimensionless arrangement. This can be prepared in a number of numerous techniques, contingent on the standards selected for distance (45) �n . F� = T ∇. n, www.nature.com/scientificreports/ and duration. These characteristics may be chosen as follows: distance and duration attributes can be chosen as R and ρ l R 3 /T , respectively. The other dimensionless quantities can be stated as follows: The procedure in considering the previous choice yields the following dimensionless numbers in the dispersion relationship:

Mathematical Formal
where − * describes the dimensionless amounts, which desired to be dropped later for a sake of easiness.
As previously stated in inequity (51) should be satisfied. Accordingly, all the following graphs are drawn in a specific domain, where the criterion (51) is automatically satisfied. Furthermore, the calculation revealed that the factor Ŵ is always positive. This indicates that the longitudinal unchanged electric field always provides a stabilizing influence, which is a preliminary conclusion. It is established by the earlier works of Melcher and Taylor 1 , Saville 2 , and many other references cited herein. At this argument, our focus is totally based on the criterion (53). Regularly, the electric field strength Log E 2 0 is often plotted versus the wave numeral k . Actually, the controlling Log E 2 0 instead of E 2 0 permits us to accumulate larger and smaller amounts at the equal vertical axis scale. It also produces smoothness of the sharp curves. In what follows, the letter S signifies the stable zones. Concurrently, the symbol U represents the unstable zones. Actually, the dimensionless procedure allows us to choose the data that has been given before, refraining from real physical problems. Consequently, the Figs. 2, 3, 4, 5, 6 and 7 depict a sample selection structure having the subsequent characteristics: Figure 2 displays the deviation of the logarithm of the dimensionless electric strength Log E 2 0 against the dimensionless wave numeral k for various values of the Laplace numeral L a . Physically, the Laplace numeral is employed to describe the free-surface fluid dynamics. It is the proportion of surface tension to momentum transfer inside a liquid. The figure shows the implication of the Laplace numeral L a . It is clear that the stable zone rises as the values L a are reduced, which demonstrates a destabilizing impact of this factor. Consequently, the Laplace numeral has a destabilizing influence on the structure. Physically, the growth in the interfacial tension, or the undisturbed cylindrical radius, or density causes more destabilization effects on the non-Newtonian cylindrical fluid. As a conclusion, for the reasonable assumptions, the non-Newtonian fluid jet is destabilized further at large amounts of the fluid Laplace numeral. Additionally, it is discovered that viscosity inhibits the emergence and growth of non-Newtonian fluid jet short-wave instability. Figure 3 demonstrates the deviation of the logarithm of the dimensionless electric strength Log E 2 0 vs the nondimensional wave numeral k of various values of the Reynolds numeral R e . As known, the Reynolds numeral     www.nature.com/scientificreports/ the non-Newtonian fluid jet. It is established that the non-Newtonian fluid jet is destabilized further clearly at great values of the Reynolds numeral for the supposed situations. Furthermore, the viscosity prevents the onset and development of short-wave stability of non-Newtonian fluid jets. Subsequently, a comparable finding was verified 50 . It should be noted that the Reynolds numeral R e helps predict flow patterns in different fluid flow situations by measuring the ratio between inertial and viscous forces. At low Reynolds numbers, flows tend to be dominated by laminar (sheet-like) flow, while at high Reynolds numbers flows tend to be turbulent. The turbulence results from differences in the fluid speed and direction, which may sometimes intersect or even move counter to the overall direction of the flow (eddy currents). These eddy currents begin to churn the flow, using up energy in the process, which for liquids increases the chances of cavitation. The Reynolds number has wide applications, ranging from liquid flow in a pipe to the passage of air over an aircraft wing. It is used to predict the transition from laminar to turbulent flow and is used in the scaling of similar but different-sized flow situations, such as between an aircraft model in a wind tunnel and the full-size version. The predictions of the onset of turbulence and the ability to calculate scaling effects can be used to help predict fluid behavior on a larger scale, such as in local or global air or A large Reynolds number indicates that viscous forces are not important at large scales of the flow. With a strong prevalence of inertial forces over viscous forces, the largest scales of fluid motion become unsuppressed and there is not enough viscosity to dissipate their motions. Since the Lorentz force is always perpendicular to the flow velocity and it impedes movement when the flow becomes turbulent (at high Reynolds numbers), the Lorentz force becomes not effective to dampen the flow, and hence the flow is destabilized as shown. Figure 4 demonstrates the deviation of the logarithm of the dimensionless electric strength Log E 2 0 vs the non-dimensional wave numeral k of various values of Weber number W e . Physically, the Weber numeral is frequently employed in the analysis of fluid flows involving two separate liquids, especially multi-phase flows through sharply rounded surfaces. It is a measure of the fluid inertia with respect to its interfacial tension. The figure can assist in understanding how a thin film flows and how droplets and bubbles improve. As a conclusion, this picture demonstrates that W e portrays a destabilizing impact on the given structure. Basically, the rise of this factor means that the primary flowing also increases. As known, the destabilizing impact of KHI is an early  www.nature.com/scientificreports/ finding confirmed by several authors cited here. Especially, this finding accords with those early obtained 50 . This depicts how amplifying the velocity, radius, or density of non-Newtonian liquid jets causes them to become unstable. Additionally, it is implied that the non-Newtonian liquid jet destabilizes more quickly with larger liquid Weber numbers. It should be noted that the Weber number W e is a dimensionless number in fluid mechanics that is often useful in analyzing fluid flows where there is an interface between two different fluids, especially for multiphase flows with strongly curved surfaces. It can be thought of as a measure of the relative importance of the fluid inertia compared to its surface tension. The quantity is useful in analyzing thin film flows and the formation of droplets and bubbles. Figure 5 indicates the deviation of the logarithm of the dimensionless electric strength Log E 2 0 against the non-dimensional wave number k of various of Darcy number. Physically, the Darcy number is important when discussing heat conduction through porous materials such as powdered metals, foams, ceramics, etc. The Darcy number in fluid dynamics due to porous media typically demonstrates the relative influence of the medium permeability versus its cross-sectional area when the diameter is squared. The Darcy number also has a considerable impact on the entropy generation caused by convection inside porous triangular additions. This chart research reveals that the unstable areas also experience. This result shows that the stable zone rises as the values D a rise, too. As shown from the mathematical formula D a , the growth in D a produces an increase in the permeability of the permeable medium. This leads to a reduction in the influence of porous media. This outcome is supported with the results of previous studies 51 . Simultaneously, the increase of D a means a decrease in the radius of the cylinder. As shown by several researchers, the increase in the cylindrical radius leads to a stabilizing influence as previously shown by Moatimid et al. 41 . Figure 6 demonstrates the deviation of the logarithm of the dimensionless electric strength Log E 2 0 against the dimensionless wave numeral k of various values of Elasticity numeral E l . Subsequently, one concludes that the elasticity numeral has a destabilizing influence on the approach that is being investigated. Physically, the increase of the relaxation duration or viscosity generates that the non-Newtonian cylindrical fluid to be more unstable. Therefore, it is concluded that, as given in the stated circumstances, the non-Newtonian fluid jet is further clearly destabilized at high values of the fluid Elasticity numeral. Furthermore, the effects of radius, or density are found to stand in contrast to the findings yielded in 52 . Figure 7 demonstrates the deviation of the logarithm of the dimensionless electric strength Log E 2 0 against the dimensionless wave numeral k of various values of the linear MHT dimensionless constant. Consequently, one concludes that α 1 has a destabilizing effect on the methodology that is being inspected. Physically, the increase in the undisturbed cylindrical radius, velocity, and MHT constant forces the non-Newtonian cylindrical fluid to be more unstable. Additionally, more heat is transferred to a point at the trough and less heat is prevented from the site as MHT increases across the interface. Therefore, it is stated in the given circumstances that the non-Newtonian fluid jet is further clearly destabilized at high values of the fluid MHT dimensionless constant. Furthermore, the interfacial tension of the non-Newtonian liquid jet is observed to avoid the beginning and improvement of short-wave stability as earlier obtained 53 .
To this end, the linear stability approach has been finished. Therefore, the next Section is depicted to display the stability configuration throughout a nonlinear approach. This approach seeks a nonlinear partial differential equation that governs the interface displacement. Fundamentally, the first step concerns with the derivation of a Ginzburg-Landau equation, which regulates the nonlinear instability standards. Consequently, the theoretical outcomes will interpret the transition curves, which classified the stability/instability zones. Apart from the linear methodology, for more accuracy, the nonlinear criteria will reveal different transition curves.

Nonlinear stability methodology
Returning to the aforementioned nonlinear characteristic equation as given in Eq. (47), the interface displacement nonlinear characteristic equation is assessed up to the third order of η to give where α(ω, k) and β(ω, k) are the second-and third-order nonlinear constants, correspondingly, and D(ω, k) is the linear coefficient.
The various time scale techniques will be used to build the analysis that follows. This methodology is used in order to demonstrate a problem-solving strategy as a function of two or additional independent parameters. The ratio of a perfect wavelength, periodic time, or modulation time scale can therefore be assumed to be a small parameter. The independent parameters z and t may be extended to include more independent parameters in the manner shown below.
while Z 0 and T 0 are both valid examples of fast variations, the slow ones are Z 1 , T 1 , Z 2 , and T 2 . The following derivative expansions can be used to express the differential operators: when θ = kZ 0 − ωT 0 , the lowest order is indicated.
The partial derivatives of time and space linear operation are (54) D(k, ω)η = α(k, ω)η 2 + β(k, ω)η 3 (55) Z n = ζ n z and T n = ζ n t n = 0, 1, 2, ........  The secular terms in Eq. (71) correspond to the coefficient of e iθ . The following solvability criterion results from the term elimination: www.nature.com/scientificreports/ It produces a complicated relationship. A uniformly valid expansion of the function η 2 takes the following form thanks to this solvability condition: By replacing both the parameters ω and k with 2ω and 2k , correspondingly, the linear dispersion function D(ω, k) can be used to structure the non-zero denominator . It must be mentioned that the disappearance of D(2ω, 2k) produces the well-known resonance curve if D(2ω, 2k) is a real function. Overall, the harmonic resonance can happen if (k, ω) and (nk, nω) meet the similar dispersion relation when n is a positive integer; for example, see Nayfeh 54 . Subsequently, the denominator transforms into a true positive function that cannot disappear. Subsequently, as will be seen later, a Ginzburg-Landau equation will be obtained.
The uniformly third order solution, which is obtained by combining Eqs. (69), (70) and (73), advances to the following solvability condition.
The second and third-order issues can then be solved in order to obtain the equation amplitude evolution. The non-security circumstances in the existence of the regularly applicable solutions in the second order, as given by Eq. (72), might be stated by using the method created 54 . According to Eq. (75), the wave approximates a second order group velocity. This indicates that the slow parameters Z 1 , T 1 through the mixture (Z 1 − V g T 1 ) affect the amplitude γ . One must arrive at the third-order equation given by Eq. in order to improve the amplitude modulation for the increasing waves (63). A single equation can be created by simplifying and combining the solvability criteria (72) and τ = ζ 2 t, The transition curves as given by the equalities of (81) and (82) should be addressed in a practical dimensionless form prior to dispersing the numerical calculations. This can be made in a variety of approaches in light of the choice of time, length, and mass appearances. Consider that the values of the factors 1/ω,a , and T/ω 2 , respectively, describe the properties of time, length, and mass.
After complex but basic calculations, the transition curve provided by Q i = 0 may be organized in a thirddegree polynomial of E 2 0 as follows.
However, the fifth-degree polynomial on E 2 0 can be used to arrange the second transition curve as given by P r Q r + P i Q i = 0 as follows: where the constants provided in Eqs. (83) and (84) are understood from the context. They will be crossed out in order to make the paper shorter.
The transition curves provided in Eqs. (83) and (84) will be investigated in order to demonstrate the stability standards used during this nonlinear stability methodology.
For this objective, consider the following particular system: In light of the above data, the calculations revealed that only one positive real root is produced by Eq. (83), however, the other two roots are complex conjugates. Really, this comes from algebra. Additionally, Eq. (84) generates only one positive real root and the other four roots in Eq. (84) are complex conjugates. Therefore, the two criteria revealed only two transition curves. Actually, the other sixth curves have no implication on the stability profile. As a surprise, the only two transition curves are coincident with each other. The symbols S 1 , S 2 indicate that the two inequalities (81) and (82) are automatically satisfied. By contrast, symbols U 1 , U 2 mean that the above stability criteria are not satisfied. Unfortunately, the nonlinear theory yields one transition curve like the linear theory. The inspection of these transition curves shows that the transition curve in the linear sense behaves like an increasing function. On the other hand, the transition curve throughout the nonlinear theory reveals a decreasing function. Therefore, one can say that the stability zone in the nonlinear sense is larger than the linear approach. Consequently, in Fig. 8, log E 2 0 is planned against the numeral of surface wavenumber k. Figure 9 demonstrates the deviation of the logarithm of the dimensionless electric strength Log E 2 0 vs the non-dimensional wave number k of various values of viscosities of liquid µ l in the nonlinear stability picture. From the examination of this picture, one concludes that the viscosity of liquid has a destabilizing impact on the approach that is being investigated. Physically, fluid resistance rises with linear viscosity and flow velocity is thus reduced, where viscosity is recognized as a factor in flow resistance. The internal fluid friction is displayed. A fluid with a high viscosity can resist motion because of the tremendous internal friction produced by its molecules. Chemically, when the intermolecular forces of attraction are strong, a material has a high viscosity. The easiest way to describe this phenomenon is to imagine two liquids rushing down a surface. This behavior is consistent with the Elasticity Number influence on the linear study previously seen in Fig. 6, where the Elasticity number is directly proportionate to the liquid viscosity. The finding is consistent with that which was previously achieved 58 .
The effects of MHT parameters α 1 and α 3 , respectively, on the stability image are shown in Figs. 10 and 11. Figure 10 demonstrates the fluctuation of the logarithm of the dimensionless electric intensity Log E 2 0 against the dimensionless wave numeral k.Consequently, α 1 has a stabilizing effect, as found in Fig. 10. The linear stability technique shown in Fig. 7 is incompatible with this influence. As seen in Fig. 11, α 3 has a stabilizing impact. Thus, it is demonstrated that the movement of heat and mass has a stabilizing effect. Additional heat may be  www.nature.com/scientificreports/ transmitted to a place where more heat is far from the position due to an increase in heat and mass transfer across the interface. As a result, it is possible to assert that the MHT factor plays a dual function in the stability standard. These findings coincide with those previously obtained by Sharma 18 . Figure 12 shows how the logarithm of the dimensionless electric field strength Log E 2 0 vs the non-dimensional wave number k with various values of permeability porosity κ varies. As shown,κ has a stabilizing effect. Physically, the porous medium permeability is a quality that quantifies the formation capacity and ability to convey fluids. Because it regulates the direction and flow rate of the reservoir fluids in the formation, the rock permeability,   www.nature.com/scientificreports/ k, is a significant rock characteristic. This effect is consistent with the impact of the Darcy numeral in the linear study as in Fig. 5, where the Darcy coefficient is directly proportional to the permeability of porosity. Figure 13 shows how the logarithm of the dimensionless electric strength Log E 2 0 against the dimensionless wave numeral k with various values of cylindrical undisturbed radius R varies. As observed,R has a dual role effect. We noted that for small values of wave numeral, the effect of radius is a destabilizing effect, whereas for larger values of wave numeral, the effect is changed to a stabilizing effect. This outcome is consistent with that which was previously achieved 8 .

Estimated surface departure formulation
The goal in this Section is to arrive at a limited solution for the surface displacement. The nonlinear technique produces the nonlinear characteristic equation that is provided in Eq. (47), as was previously demonstrated. With the complex constants of the interface deviation η(z, t) , it forms a nonlinear second-order differential equation. Actually, this equation analysis in its current form is quite challenging. Physically, the amplitude departure η(z, t) nature must be a real function. One may only consider the time dependent, i.e.,η = η(0, t) = Ŵ(t) , in order to make the subsequent calculations easier, the assumption is that It will be easier to perform the calculations that follow if only the time dependence ( z = 0 ), which is the true indicator of historical stability, is considered. It is possible to rewrite the nonlinear characteristic equation Eq. (47) as w h e r e G(η) = i f 0 η + (g 1 +i f 1 ) η t + (g 2 +i f 2 ) η 2 + (g 3 +i f 3 ) η 3 + (g 4 +i f 4 ) η η t +(g 5 +i f 5 ) η 2 η t +g 6 ηη tt + g 7 η 2 η tt .
(86) η tt (t) + (a 1 + i (b 1 + ka 4 )) η t + (a 3 − kb 2 − k 2 a 5 ) + i (ka 2 + b 3 ) η + G(η) = 0,  where the factors ϑ j will be calculated later on after being combined with the initial features of the issue at hand. Disregarding the secular terms, this will be accomplished. It follows that the Homotopy equation might be formulated as follows: where the symbol ω 2 = (a 3 − kb 2 − k 2 a 5 ) is the present model natural frequency. This method is primarily dependent on a small factor and is based on an expanded frequency analysis. This method allows the enlarged artificial frequency ̟ 2 to be stated as follows 59 : The expanded form of the time-dependent function η(0, t, ) is as follows: By combining Eqs. (89) and (88), the nonlinear characteristic equation can be rewritten as Considering the Laplace transform applied to Eq. (91) and the initial conditions specified in Eq. (87), the result is: Employing the inverse Laplace transforms of both sides of Eq. (92), one finds Utilizing Eq. (90) improvement the dependent parameter η(0, t, ) and then locating the constants with equal powers on both sides, one is able to obtain We will substitute from Eq. (94) into Eq. (95) and remove the secular terms. For this, it is better to ignore the coefficients of the functions cos ̟ t and sin ̟ t at all phases; solving the two equations simultaneously produces the values of R 1 and ϑ 0 .
The periodical solution at this point in time is provided by (87) η(0, 0) = 0, and η t (0, 0) = ∞ j=0 j ϑ j , : In reality, the restricted estimated solution in Eq. (100) requires that the trigonometric function arguments be real numbers. It can be concluded that the extended frequency fulfills a particular characteristic equation when Eqs. (89) and (97) are combined. The analysis revealed that, at the extended frequencies, this equation resembles a polynomial of the second degree. This equation can be solved to determine the value.
The following figure is a graph of a structure receiving the following specifics: The computations revealed that the value of the expanded frequencies exists and equals ̟ = 2.69 . The disturbing response from Eq. (100) will be plotted in Fig. 14. This figure demonstrates the periodic behavior of the profile of the surface displacement. The prior behavior results from the secular terms of the non-homogeneous parts of the controlling equation of motion of the surface waves that is cancelled out.

Concluding remarks
In this study, we have examined the linear and nonlinear stability of a vertical cylindrical interface that separates between two homogeneous incompressible dielectric viscoelastic liquids. The inner and outer liquids are occupied by a non-Newtonian Oldroyd-B type. The Oldroyd-B fluid possess both relaxation, retardation times in addition a viscosity parameter. The system is acted upon by an unchanged longitudinal electric field. A basic construction of VPF is utilized to reduce the mathematical workload. Numerous specific cases are represented based the selected data. Through porous media, Darcy's law controls the flow. The permeability of the two fluids is referred to simply as for convenience. Darcy's law, an experimental hypothesis that describes the creeping flow of Newtonian liquids through porous media, is supported by scientific evidence. Additionally, MHT is reflected with the shortened formulation of Hsieh 38 and 39 . The nonlinear approach contributes primarily to the solution of the fundamental linear equations of motion and using the relevant suitable nonlinear border circumstances. This process produces a nonlinear characteristic equation of the surface displacement of the surface waves. The dimensionless procedure exposes a collection of dimensionless physical numbers, in addition to our prior work, see Moatimid and Mostapha 60 . A linear dispersion relation is obtained by disregarding the nonlinear elements. The linear stability standards are judged via the Routh-Hurwitz theory. To confirm the theoretical outcomes, numerical estimations are carried out. A series of diagrams are depicted in order to confirm similar findings to those reported in the literatures. For the nonlinear characteristic equation, a novel approach is utilized to achieve the stability benchmarks. Consequently, the stability standards are theoretically and numerically validated. Our upcoming study will address these forms non-Newtonian models in consideration due to the practical applications of the various structures of non-Newtonian prototyping as progress works. Furthermore, particular linear (99) ϑ 1 = ϑ 2 0 96̟ 6 f 1 + 72̟ 4 g 5 ϑ 2 0     32̟ 4 f 1 g 4 + (9̟ 2 f 3 g 1 − 9̟ 2 f 1 g 3 + 24̟ 2 g 2 g 4 + 9̟ 4 f 1 g 7 )ϑ 0 + (48f 2 f 3 − 48g 2 g 3 +24̟ 2 g 4 g 5 + 48̟ 2 g 3 g 6 + 48̟ 2 g 2 g 7 − 48̟ 4 g 6 g 7 )ϑ 2 0 + (12f 3 f 5 − 12g 3 g 5 + 18̟ 2 g 5 g 7 )ϑ 3 0 +̟ 2 f 0 (64f 2 + 3f 5 ϑ 0 ) − 8̟ 2 f 4 (4̟ 2 g 1 + 3(ϑ 0 f 2 + f 5 ϑ 2 0 )) + (64̟ 2 g 2 − 64̟ 4 g 6 + 3̟ 2 g 5 ϑ 0 )ℜ 1  www.nature.com/scientificreports/ stability non-Newtonian standards were explored along with our earlier articles 10,46,[50][51][52] . Therefore, away from the simplified formulation of Hsieh's formula 38 and[39[, the future works will adopt the energy as well as the concentration equations. Without any doubt, this approach reveals several parameters concerning the MHT phenomenon. The nonlinear stability approaches will be emphasized in order to get a higher degree of precision. Overall, the stability zones are enhanced in the nonlinear instability approach in distinction with those in the linear approach. For more convenience, the stabilizing/destabilizing influences of various physical parameters in the problem at hand, may be summarized throughout the following: • • Additionally, throughout the nonlinear stability methodology, the graph of the surface displacement has been established in a uniform formula after the elimination of the secular terms.
As we make progress, the following topics will be addressed in later papers: • In addition to providing the existence of the energy and concentration equations, Hsieh's simple approach will no longer be used. • Various nanofluids with the numerous non-Newtonian fluids will be implemented.
• Different stability issues in a plane geometry will be addressed in light of the importance of numerous actual engineering implementations. • A time-varying exterior fields will be supplied.
• An innovative method will be used to examine the stability methodology in light of the non-perturbative approaches in evaluating the nonlinear partial differential equations..